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ABSTRACT 

In this paper we design and analyze algorithms for asynchronously 
solving linear programs using nonlinear signal processing structures. 
In particular, we discuss a general procedure for generating these 
structures such that a fixed-point of the structure is within a change 
of basis the minimizer of an associated linear program. We discuss 
methods for organizing the computation into distributed implemen¬ 
tations and provide a treatment of convergence. The presented algo¬ 
rithms are accompanied by numerical simulations of the Chebyshev 
center and basis pursuit problems. 

Index Terms — asynchronous optimization, distributed opti¬ 
mization, linear optimization, nonlinear signal processing 

1. INTRODUCTION 

Traditional linear programming algorithms typically iterate an or¬ 
dered sequence of operations until either a suitable stopping crite¬ 
rion is met or a minimizer is identified, e.g. interior point and basis 
exchange methods [1]. In transforming such an algorithm into the 
distributed setting, issues related to task allocation, communication, 
and global versus local synchronization become of paramount im¬ 
portance. Computational reorganization often involves partitioning 
the sequence of operations into those which distribute in the sense 
that their execution can make use of many processors that do not 
interchange information and those which do not, e.g. methods like 
[2]-[4] and the distributed simplex method in [5]. Global optimiza¬ 
tion algorithms in the style of monte carlo methods paired with lo¬ 
cal search techniques make efficient use of distributed resources but 
often require careful parameter tuning for a given problem [ 6 ]. Un¬ 
derstanding the complexity associated with a sequential algorithm 
differs in many ways from its distributed counterpart as the penalties 
incurred for communication and synchronization may be more ex¬ 
pensive than that suggested by traditional complexity measures [7]. 

In this paper we specifically address linear programs from a con¬ 
servative signal processing perspective consistent with the general 
framework presented in [ 8 ]. In particular, we design and analyze a 
generally nonlinear signal-flow structure which when implemented 
asynchronously solves an associated system of stationarity condi¬ 
tions that was developed in [9]. The viewpoint of solving either 
an optimization or constraint satisfaction problem using an asyn¬ 
chronous signal processing system naturally lends itself to under¬ 
standing important issues such as algorithmic scalability, robustness 
with respect to communication and processing delays, and compu¬ 
tational heterogeneity. Furthermore, issues pertaining to the identi¬ 
fication of sufficient conditions for convergence may be addressed 
using well-known stability and dynamic systems results. 
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2. PRELIMINARIES 

The optimization problem of minimizing a linear objective function 
subject to linear equality and inequality constraints is described in 
standard form as 

minimize 
^ ~ 
subject to Ax < b 
£ > 0 

where x £ is the decision vector, / € is the cost vector, 
b G R“^ is the constraint vector, and A G R“^^^ is the coeffi¬ 
cient matrix. We describe an asynchronous algorithm for solving 
( 1 ) in this paper via a signal-flow structure consisting of subsys¬ 
tems realized as maps coupled together using asynchronous delays, 
i.e. randomly triggered sample-and-hold elements which we model 
as independent Bernoulli processes. The general strategy underly¬ 
ing the presented algorithms described in this way is to determine a 
solution to a system of stationarity conditions associated with a par¬ 
ticular reformulation of ( 1 ) by interconnecting linear and nonlinear 
signal-flow elements such that they collectively describe the behav¬ 
ior of the stationarity conditions, resolving any delay-free loops, and 
finally running the system to a fixed-point. The delay-free loops are 
either resolved algebraically, e.g. using the automated techniques in 
[ 10 ], or by inserting asynchronous delays depending on the specific 
form of the associated nonlinearity. 

We refer to a linear program problem statement organized ac¬ 
cording to the following conventions as being in asynchronous form: 

(i) every vector (except the cost vector) is involved in a system 
of linear equations, 

(ii) every vector in (i) is either fixed, unconstrained, or non¬ 
negative, 

(iii) every vector in (i) with non-zero cost coefficients must be un¬ 
constrained. 

Indeed, a linear program in standard form can always be recast into 
asynchronous form as 

minimize 

subject to Bzj = 23 ( 2 ) 

^2 >0 

where the minimization is explicitly over those variables which are 
unconstrained and non-negative and 
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Note that we have introduced a non-negative vector y G in 
order to enforce the linear inequality Ax < b and have made two 
equality-constrained copies of x where Xj encodes the non-zero cost 
coefficients and Xj enforces the non-negativity constraints. 











Consistent with the presentation in [8][11], the system of sta- 
tionarity conditions for the formulation in (2) is of the general form 

d* = Gc* (3) 

c* = rn{d*) (4) 

where G is an orthogonal matrix, m(-) is an element-wise memory¬ 
less nonlinearity, c and d denote respectively the input to and output 
from G, and c* and d* denote a fixed-point of the algebraic system. 
The vectors c and d are in particular organized according to the or¬ 
dering of the inputs followed by the outputs of the linear equality 
constraints in (2), i.e. 
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Furthermore, the matrix G in (3) is generated using B in (2) as 
G = iI + R){I-Ry^ 


(5) 
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where R — —R^ is a skew-symmetric matrix of the form 
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Fig. 1. A general signal-flow structure (top left) used in the initial 
description of the stationarity conditions for linear programs in asyn¬ 
chronous form and the resulting signal-flow structure (top right) in¬ 
dicating three possible locations for distributing system state specif¬ 
ically indicated using dashed boxes. 


in asynchronous form consistent with the general strategy previ¬ 
ously described. We begin this presentation by partitioning the lin¬ 
ear stationarity conditions in (3) where we specifically delineate be¬ 
tween those variables associated with memoryless nonlinearities in 
(4) which are affine maps and those which are generally not, i.e. 


The orthogonality of G is readily verified using the fact that I + R 
and 7 — 77 commute. Moreover, it follows that G is an element of the 
subset of special orthogonal matrices which do not have eigenvalues 
of —1. When B is itself an orthogonal matrix, e.g. the discrete 
Fourier transform matrix, the expression in (5) simplifies to G = 77 
and thus if a fast implementation of B is available it may be used in 
the implementation of G. 

The memoryless nonlinearities used to realize the stationarity 
conditions in (4) for a linear program described in asynchronous 
form are listed in Table 1 for a single element « which either be¬ 
longs to (input) or (output). The nonlinearities are specifically 
described for the various pairings of cost contributions and set mem¬ 
berships which arise in linear programs cast this way. For example, 2 
being unconstrained with no contribution to the overall cost follows 
from the second row with p = 0. In the sequel we resolve delay-free 
loops in the presented signal-flow structures associated with the first 
two rows algebraically and the third row using asynchronous delays. 

Table 1. nonlinearities for linear programs 


cost 

set membership 

m(-) (input) 

m(-) (output) 

none 

fixed z = p 

c = —d + 2p 

c = d — 2p 

pz 

unconstrained 

c = d — 2p 

c = —d + 2p 

none 

non-negative 

c= |d| 

c = -\d\ 


Given a fixed-point c* and d* of the stationarity conditions (3)- 
(4), the argument a;* which minimizes (2) and consequntly (1) may 
be determined by extracting either a; or from 

it = i(4+4) ^ ^2 = 5(4-4)- (6) 

3 . ASYNCHRONOUS LINEAR PROGRAM ALGORITHMS 

In this section we present a general asynchronous algorithm in the 
form of a signal-flow structure for solving linear programs described 
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where the variables are ordered according to 
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and where G is generated using (5) and block partitioned accord¬ 
ingly. Eliminating those variables associated with b and in (3), 
i.e. dj and , results in an affine system of the form 

^2 ~ ^ T2 7 ” — 

where 

G' = G22-FG2i(7-S'Gii)"^SGi2 

e = 2 G 21 (T-BGii)-^ 

Q _ —7m 0 

I. 0 7 iv • 

The stationarity conditions corresponding to the described reduced 
representation are given by (8) paired with y = “ 1 ^ 21 - ^ signal- 
flow structure depicting a fixed-point of these stationarity conditions 
is portrayed in Figure 1 on the top left. Another class of related 
algorithms follow from defining a sequence of signal-flow systems 
which smoothly deform into this structure. Referring again to Fig. 1, 
the signal-flow structure on the top right depicts three possible loca¬ 
tions at which asynchronous delays may be inserted in order to break 
delay-free loops where the dashed boxes labeled D denote a vector 
asynchronous delay element. This structure forms the basis from 
which various implementations may be synthesized. Seven natural 
organizations of the system state are additionally depicted. 
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Signal-flow structure storage 
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Fig. 2. Left: a signal processing system organized into a conceptual 
associative array. Right: the general computation for the fc* update. 


4. DISTRIBUTED ASYNCHRONOUS IMPLEMENTATIONS 

We next briefly present an example implementation of an asyn¬ 
chronous signal-flow structure using a conceptual associative array 
organized using the key-value pairs depicted on the left of Fig. 2. 
The algorithms used in Section 6 were specifically implemented us¬ 
ing this approach where the general computation associated with an 
asynchronous update is depicted on the right and explained below. 

Consider the causal system, i.e. with initial condition CjlO] = 0, 
defined by the recurrence relation 

C2N = rn{d^[n-l]) 

dj [n] = G'cj [n] + e ^ ^ 

and note that a fixed-point of this system is a fixed-point of (8). For 
n > 1 an equivalent description obtained via manipulation is 

K 

c*2N = (»^(rf2.fch- 1]) - C2,fc[«'- 1]) (13) 

k^l 

where denotes the fc* column of G' and (•)2,fc denotes the fc* 
element of (•)2. The initial condition ^2(0] = e is required for (12) 
and (13) to produce the same output d2[n] for all n > 0 . An asyn¬ 
chronous implementation of this system then follows from compu¬ 
tational nodes executing the general computation described in Fig. 2 
on the right for a randomly selected k. The “compute” stage com¬ 
putes the fc* term of the summation in (13) using information ob¬ 
tained in the “lookup” stage, while the “increment” stage is used to 
update dj and Sa k without requiring a full read-and-write operation. 

5. CONVERGENCE ANALYSIS 

In this section we analyze convergence properties of the signal-flow 
structure in Fig. 3 on the left for both synchronous and asynchronous 
settings. To facilitate this analysis, let T denote an iterated operator 
mapping the output of the delay to the input such that composing T 
with itself ad infinitum produces a synchronous implementation, i.e. 

^(^2) ~ G'm{d 2 ) + 6- (14) 

Indeed, for the given operator linear convergence to a unique fixed- 
point is guaranteed provided that T is Lipschitz continuous with con¬ 
stant Lt G [0, 1), i.e. for every d® and dj^^ T satisfies 

• (15) 

Establishing the existance of a fixed-point for this case follows im¬ 
mediately by assigning dj^^ = d2[n] and d^^^ = d2[n — 1] and 
taking the limit of the iterated inequality, i.e. 

lim ||d2 [n] - d2 [n - IJHj < ||4 [1] - d2 [0]||2 lim = 0. 

n—Foo n—Foo 

Suppose two signals dj^l [n] 7^ dj^l [n] are both fixed-points of T. 
Taking the limit of the iterated inequality in (15) for this case yields 

lim lldj^^ [n] — d^^ [n]|| < lld^^^ [0] — d^'' [0]|| lim LJ = 0. 
n—Foo II II2II ll2'^~^c>o 



Fig. 3. The signal-flow structure used for convergence analysis. 

Since a fixed-point is defined as being independent of n we have 
dj^l [n] = d^^ [n] which contradicts the original assumption, thereby 
establishing that the fixed-point of T is unique. 

We proceed to analyze the asynchronous setting using the sys¬ 
tem in Fig. 3 on the right consisting of the difference between the 
system on the left and a fixed-point. Note that when m(-) satisfies 
(16) with Lipschitz constant Lm then the map m' from d^ [n — 1] to 
Cjfn] does too. Let An-m denote the event that a particular asyn¬ 
chronous delay element last fired at time n — m and let d[m] denote 
the corresponding geometric distribution with parameter p € (0,1). 
Then, using the law of total expectation on Cj [n] yields 

00 

m=0 

00 

m=0 

Summing the K elements composing C 2 , interchanging the summa¬ 
tions, and rewriting the expression using vectors yields 


00 

[ll&Nla] = (4[n-m-1])||2](16) 

m =0 

00 

< Em ^ - ”1-- IJlIa] (IE) 

m =0 

where we have explicitly used the Lipschitz continuity of m! in go¬ 
ing from (16) to (17). Finally, rearranging indices and noting that 
E [lld^Wlli] = E [||c 2 [n]|||i leads to the expression 
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which takes the form of a convolution bound and is sufficient for 
asymptotic average convergence if Lm € [0,1). 

For linear programs, however, T is readily verified to satisfy 
Lt = 1 since G' is orthogonal and m(-) is non-expansive. Justify¬ 
ing the orthogonality of G' in (9) follows from showing that G' is a 
linear isometry. In particular, since G and S are orthogonal matrices 


- rfi 


( 2 ) 


\\di^^ -c *2 


( 2 ) 




2 

2 


and since c, = Sd, + k for some constant k it follows that 





Since the operator T as described is on the boundary of guaranteed 
convergence, it stands to reason that a homotopic relaxation may 
remedy divergent or oscillatory behavior. Toward this end, define an 
operator T' with homotopy parameter a € [0, 1] as 


T'{d 2 , a) = aT{d 2 ) + (1 - a)To{d 2 ) 


where To is chosen such that Ltq G [0, 1) and varying a from 0 to 
1 corresponds to a smooth deformation from To to T. Therefore, it 
follows that synchronous or asynchronous implementations of T' as 
a —> 1 converge to a fixed-point. 
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Convergence for the Chebyshev center problem 
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Fig. 4. Asynchronous convergencefor the Chebyshev center problem 
averaged over 500 trials. 


6. NUMERICAL SIMULATIONS 


In this section the convergence properties associated with asyn¬ 
chronously solving two linear programming problems using the 
algorithms developed in this paper are explored where the asyn¬ 
chronous delays were numerically simulated using discrete-time 
sample-and-hold elements triggered by independent Bernoulli pro¬ 
cesses. For the sake of comparison, we use the metric of an equiva¬ 
lent iteration to normalize between various probabilities of sampling, 
i.e. an equivalent iteration is the same total amount of computation 
associated with a synchronous iteration where all delays fire inde¬ 
pendent of the probability of sampling. 


6.1. The Chebyshev center problem 


Consider as an example the Chebyshev center problem [12] given by 
minimize —r 

i£c ’ ^ 

subject to + llfflilhr <bi 1 < i < M 

r > 0 

In this form (18) is explicitly identifying the largest Euclidean ball 
which can be inscribed within a convex polytope described in half¬ 
space representation by {z: Az < b} where r and denote the 
balls radius and center, respectively, and a^ is the i* column of A. 
We proceed by recasting this problem into asynchronous form as 


minimize 

■£c’— 

subject to 


-n 


■ 1 0 o' 


ri 


r2 

n —A b 
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r2 > 0, 2 > 0, t = 1 


(19) 


where the i* entry of n is ||ai||2. The organization of variables for 
this problem along with the memoryless nonlinearities associated 
with the stationarity conditions in (4) is summarized in Table 2. The 
matrix in (19) is used to generate G using (5) and in the process of 
synthesizing an asynchronous algorithm adhering to the presented 
framework the system variables associated with ri, and t are 
eliminated using (8) for an appropriate choice of dj, Cj, dj- snd 1-2• 
Continuing with this notation, the balls center x* may be recovered 
from cj by first recovering Cj and dj; via 

dt = (7-Gii5')-^(Gi2C*-FGiib) (20) 

Cj = Sdl + h (21) 


followed by partitioning z* in (6) where 


h = 


2 

0 

2 


and S = 


I 0 
0 -1 


( 22 ) 



Fig. 5. Asynchronous convergence for the basis pursuit problem 
averaged over 1000 trials. 

Convergence results for the Chebyshev center problem in (19) 
averaged over 500 runs of an asynchronous algorithm are depicted 
in Fig. 4 on the right where the convex polygon is defined using 
200 random hyperplanes in a 100-dimensional space. As a further 
example, the geometric figure on the left tracks the center of the 
Euclidean sphere in 3-dimensions for the depicted polygon over the 
course of an asynchronous algorithm as it converges to its fixed- 
point. The center is tracked using a blue line and the final Euclidean 
sphere is also depicted. 


Table 2. variable organization for (19) 


variable 

cost 

set membership 

I/O 

implementation 

ri 

—n 

unconstrained 

input 

Cri — 2 


none 

fixed 

input 


t 

none 

fixed 

input 

Ct = —dt + 2 

r2 

none 

non-negative 

output 

Cr2 — “1^7-2 1 

£ 

none 

non-negative 

output 



6.2. The basis pursuit problem 


Consider as another example the basis pursuit problem given by 


minimize ||®||r 
subject to A® = b 


(23) 


Although the process of recasting this problem as a linear program 
in standard form is well known, we proceed by using a nonlinearity 
which corresponds to an unconstrained variable z with cost contri¬ 
bution \z\. In particular, for « being an input to the problems linear 
equality constraints the associated nonlinearity is 

\ d+ 2, d<-l 

c = mi(d) = < —d, |d| < 1 

[ d-2, d> 1 

and for z being an output the nonlinearity is c = —mi{d). Conver¬ 
gence results for a 16-sparse x* in a 512-dimensional space recov¬ 
ered using 200 random measurements averaged over 1000 runs of 
an asynchronous algorithm adhering to the presented framework are 
depicted in Fig. 5. In particular, the objective value, log of the dis¬ 
tance between d^ to d*, and log of the distance between x to x* are 
depicted as a function of equivalent iteration where the delays fire 
with probability p = 0.2, 0.4, 0.6, and 0.8. A homotopic method 

was used to encourage initial convergence, i.e. the nonlinearity used 

,2 

for equivalent iterations k — 1,..., 10 was c = (1 — 0.95 )mi (d) 
and c = mi (d) thereafter. 





































































7. REFERENCES 


[1] D. Bertsimas and J. N. Tsitsiklis, Introduction to Linear Opti¬ 
mization, Athena Scientific, 1997. 

[2] N. Parikh and S. Boyd, “Block splitting for distributed opti¬ 
mization,” Mathematical Programming Computation, 2014. 

[3] P. A. Forero, A. Cano, and G. B. Giannakis, “Consensus-based 
distributed support vector machines,” J. Mach. Learn. Res., 
2010 . 

[4] E. Wei and A. Ozdaglar, “Distributed alternating direction 
method of multipliers,” in Decision and Control (CDC), 2012 
IEEE 51st Annual Conference on, Dec 2012, pp. 5445-5450. 

[5] M. Burger, G. Notarstefano, F. Allgower, and F. Bullo, “A 
distributed simplex algorithm and the multi-agent assignment 
problem,” m American Control Conference (ACC), 2011, June 
2011, pp. 2639-2644. 

[6] T. Desell, Asynchronous Global Optimization for Massive- 
Scale Computing, Ph.D. thesis, Rensselaer Polytechnic Insti¬ 
tute, 2009. 

[7] E. Solomonik, E. Carson, N. Knight, and J. Demmel, “Trade¬ 
offs between synchronization, communication, and work in 
parallel linear algebra computations,” Tech. Rep. UCB/EECS- 
2014-8, EECS Department, University of California, Berkeley, 
Jan 2014. 

[8] T. A. Baran and T. A. Lahlou, “Conservative signal processing 
architectures for asynchronous, distributed optimization part 1: 
General framework,” in Proc. of IEEE Global Conference on 
Signal and Information Processing (GlobalSIP), 2014. 

[9] T. A. Baran, Conservation in Signal Processing Systems, Ph.D. 
thesis, Massachusetts Institute of Technology, 2012. 

[10] T. A. Baran and T. A. Lahlou, “Implementation of intercon¬ 
nective systems,” in Proc. of IEEE International Conference 
on Acoustics, Speech, and Signal Processing (ICASSP), 2015. 

[11] T. A. Baran and T. A. Lahlou, “Conservative signal processing 
architectures for asynchronous, distributed optimization part 
11: Example systems,” in Proc. of IEEE Global Conference 
on Signal and Information Processing (GlobalSIP), 2014. 

[12] S. Boyd and L. Vandenberghe, Convex Optimization, Cam¬ 
bridge University Press, New York, NY, USA, 2004. 



